The spectral form factor is not self-averaging 
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The spectral form factor, k(t), is the Fourier transform of 
the two level correlation function C(x), which is the averaged 
probability for finding two energy levels spaced x mean level 
spacings apart. The average is over a piece of the spectrum 
of width W in the neighborhood of energy eo- An additional 
ensemble average is traditionally carried out, as in random 
matrix theory. Recently a theoretical calculation of k(t) for 
a particular system, with an energy average only, found in- 
teresting nonuniversal semiclassical effects at times t ~ 1 in 
units of r% = 2irh/(mean level spacing). This is of great in- 
terest ifk(t) is self-averaging, i.e, if the properties of a typical 
member of the ensemble are the same as the ensemble average 
properties. We here argue that this is not always the case, 
and that for many important systems an ensemble average is 
essential to see detailed properties of k(t). In other systems, 
notably the Riemann zeta function, it is likely possible to see 
the properties by an analysis of the spectrum. 



Recent work connects the energy level statistics of 
a particular system, say a Bunimovich stadium billiard, 
with the predictions of random matrix theory [RMT] [^j . 
This work has generated great interest and some contro- 
versy. In particular it was suggested that the form factor 
k(t) of the spectrum of a given system can be calculated 
theoretically to good approximation directly in terms of 
certain classical parameters [i.e. the periodic orbits] of 
the system. For scaling systems like billiards, not only 
does the form factor approach the universal RMT result 
in an appropriate high energy limit, but predictions [0, 
H are made as to how this limit is approached. These 
predictions differ in detail from one another, but their 
gross features are the same. 

It would obviously be of interest to take the spectrum 
of some system, obtained either numerically or experi- 
mentally, and use it to calculate k(t), so that comparison 
of the approximate theories with the numerically 'exact' 
structure factor can be made. The main result of this 
paper is that in many cases such a calculation is unex- 
pectedly difficult, and may even be impossible in principle. 
We refer to this difficulty as the fact that the form factor 
is not a self- averaging quantity. 

Hard chaotic billiards without time reversal symme- 
try, 'GUE' ||, the case worked out in detail up to now, 
unfortunately fall into the probably impossible category. 
Billiards with nonisolated orbits are likely to be difficult 
although not impossible to study numerically. However, 
the Riemann zeta function, on the hypothesis that its ze- 



roes are eigenenergies of some linear Hermitean operator, 
can probably be studied effectively, especially in view of 
the large database of zeroes obtained by Odlysko Q| . 

Of all the two level correlation functions, [which are 
interrelated by linear transformations], the form factor 
is of particular interest because, according to the above 
ideas [which originated with Berry Q], the dimension- 
less time t, the argument of k(t), is essentially the rel- 
evant classical orbit period in units of the Heisenberg 
time, t% = 2ir%/ A, where A is the mean level spacing. 
In particular, for small times t, it has long been estab- 
lished [and exploited ||] that the form factor is strongly 
peaked at the periods of the short periodic orbits, and 
the weights of these peaks are determined by the stabil- 
ity properties of the orbits. Other correlators smear over 
the time and are not so directly interpreted classically. 

The new feature of the above theories is that there are 
definite nonuniversal effects in a particular range of t, 
namely near t — 1. These effects disappear in the high 
energy limit, as some power of energy. Exactly which 
power determines whether the effects can be observed in 
numerical calculations. 

The form factor is usually loosely defined as the Fourier 
transform of the two level correlation function C (x) , 

C{x) = A 2 H^5{E + \xA - E a )S(E -±xA-E b )\ 

\ a,b I w 

(1) 

where the spectrum is given by the sequence of levels 
E a . The average of E is over an energy range W about 
a central energy eo. An additional average historically 
has been over an ensemble of systems, for example ones 
given by RMT or by a random scattering potential. We 
are investigating the consequences of disposing of this 
second average. The form factor is then 



k{t) = dx 



(2) 



Two quite different techniques have been used. One 
|l| starts from a supersymmetric integral representation 
of C{x), and then makes some approximations which are 
discussed and justified. The other || extends the diago- 
nal approximation |5| in an imaginative way to all values 
of t. In both cases it is assumed that the energy average 
by itself validates the operation used to give the diag- 
onal approximation, [DA], [for t << 1] which is equiva- 
lent to the idea that the actions of all orbits [unrelated 
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by symmetry] may be separately averaged. In addition, 
statistics of classical orbits are invoked in the form of the 
Hannay-Ozorio de Almeida sum rule , and possible ex- 
tensions of it (!]] . 

A third technique [p| explicitly combines random ma- 
trix theory with orbital information about a particu- 
lar system. It therefore does make an ensemble average, 
something like an average over the random class of bil- 
liards all with the same area and also whose short orbits 
have the same classical periods and stabilities. The re- 
sults are less complete, but do agree qualitatively with 
the other methods. 

Ensemble averages are often used as a calculational 
device, even when no ensemble is actually contemplated. 
The resistance of a particular copper wire is conveniently 
calculated by averaging over an ensemble representing 
the 'random' position of impurities. An experimentalist 
can check the theory by measuring the resistance of just 
one wire, however. That's because resistance properties 
of a typical member of the ensemble are the same as the 
average resistance over the ensemble: resistance is in this 
case a self- averaging quantity. Self-averaging is similar 
to ergodicity. For an ergodic system, a time-averaged 
quantity is self-averaging. 

Clearly self-averaging is a valuable trait. It is therefore 
disappointing to find that the spectral form factor, by one 
standard definition, is not self-averaging. Other defini- 
tions can be used, and some of them are self-averaging. 
However, the alternative definitions known to us are de- 
fective in the sense that in important cases they eliminate 
the predicted effects: they are not sufficiently flexible to 
extract the small predicted effects from the spectral data. 

The energy average is essential, [even with an ensemble 
average] , but we are interested in the case that the results 
depend only weakly on the averaging window W , and 
indeed, they should not depend much on exactly how the 
average is made. This is supposed to be the case if the 
following two facts hold. 

• There are a large number of energy levels in the 
average, i.e. w = W/A >> 1. 

• The averaging width W is classically small. This 
means that energy dependent objects like the mean 
level spacing, and other quantities characterizing 
the underlying classical system, do not vary much, 
[or can be linearized] , in the energy range eo ± h W. 
For short we describe this as W << eo, which tac- 
itly supposes a natural zero of the energy. 

Theorists like to clean things up by taking a limit: 
e /W — ► oo, W/A — > oo. This requires a system with an 
infinite number of levels, and also requires knowledge of 
how things change as eo gets very large. A prominent 
example of such a system is that of billiards. The order 
of limits is not completely trivial however and another 
parameter is needed. Such limits should be left to the 
very end, and then scrutinized carefully. In fact, numeri- 



cal calculations can't take such limits, so one must study 
how they deviate from the limiting behavior. 

The energy average alone certainly suffices to give in- 
teresting results for small t, t ~ t c , where t c is the scale 
of the periods of the shortest periodic orbits measured in 
units of the Planck time. For such short times, k (t) gives 
the length spectrum [^, that is, it has peaks at t = T p /th, 
i.e. at the periods of the short periodic orbits, which for 
billiards are proportional to the length of the orbits. 

Here is a more careful definition. Let the energy 
average be given by a Gaussian, i.e. (f(E)) w = 
J dEf{E)exp( : -4:{E~e () ) 2 /W 2 ) . We argue next 

that the x integral in Eq.(||) should be cut off, and we 
choose a Gaussian: exp(— (xA/W x ) 2 ). We claim that a 
natural choice is W x ~ W, and the particular choice 
W x = W is especially congenial. 

Clearly W x should not be too large, since approxima- 
tions relying on xA << eo must be made. If W is as 
large as possible, then W x can't be much bigger. There 
is no advantage to taking W unnecessarily small, as one 
would be throwing away data, in a numerical calculation, 
for no gain. So, we take W x < W. 

The choice W x « W has virtue. However, the same 
effect can be achieved by making a time-convolution on 
k(t), i.e. smearing each t by an amount r. We find it 
convenient to choose this route. Therefore, we define 



k(t) 
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where C is based on the Gaussian energy average. 

This k{t) [except for |t| << 1] does not look like the 
ensemble average. Rather, it has very considerable 'noise' 
around the ensemble mean value. 

To see this, we remark that k(t) can be written 



k(t) 



F a e 



(4) 



where F 2 = exp(-4(.E a - e ) 2 /W 2 ). The partic- 

ular choice of cutoff W x = W has the advantage that 
the absolute square of a single sum appears, rather than 
requiring a double sum. It indicates, for example, that 
k(t) should be close to a positive function, no matter the 
form of the averaging. More general functions, F ai ap- 
propriately normalized, will do just as well. Note that for 
large t, one expects just the diagonal terms in the double 
sum to survive, so that in some sense k — ► 1 in this limit. 
At very small t, it is expected that k(t) ~ 5(t). But in 
this case, the sum can be replaced by an integral, and a 
Gaussian of width A/W replaces the 5 function. 

Eq. (|J) makes it clear that for the spectrum of a chaotic 
system, k(t), for large t, is approximately given by a ran- 
dom walk in the complex plane, of about w steps. For 
the kind of systems envisaged, there will be significant 
random walk character even for t ~ 1. Insofar as the 
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phases of the individual terms in the sum Eq. (J4|) are in- 
deed independent random variables, one may show that 
the distribution of values of y — k(t) is 

p(y) = e- y . (5) 

This result is independent of the form of F a and of W, as 
long as there are many terms in the sums of Eq.(Q), i.e. 
as long as w >> 1. Basically, the result is an application 
of the central limit theorem. The mean value [over an 
ensemble of such random variables] is (k(t)) = 1, and the 

variance of the distribution is (^(k(t) — l) 2 ^} = 1. 

For smaller t, think of the E a 's as ordered. Then it will 
take some number of steps, of order g(t), say, before the 
knowledge of the original phase is lost. There will then 
be w/g(t) effective random walk steps. The effective step 
length is still of order unity, because the phase changes 
systematically by an amount e 2lxlt at each step. This 
gives a formula for the distribution of k(t), p(k(t)) = 
g(t)e~ k ^ 9 ^ . A natural guess for g(t) gives 

k E {t) 

where ksit) is the ensemble average form factor. 

Thus, for any averaging scale W, the calculated value 
of k(t) will suffer large fluctuations, when t is of order 
unity or greater. Increasing W does not change the dis- 
tribution, but it does make k(t) vary more rapidly. It's 
clear that if t is changed by an amount of order A/W then 
k(t) changes appreciably. For example, one estimates the 

random walk average of ^ | ^ J2 a F a e 27rit ^ Ea -~ e °' )/A = 

From this it follows that smearing the time t by an 
amount r is like averaging over a number tW/A inde- 
pendent choices from the distribution p(k(t)). Thus the 
'noise' in the smeared function k T (t) is reduced by a fac- 
tor y/lJwT and of course k T (t) changes appreciably only 
when t is varied by amount of order r. 

There is literature which numerically displays the noise 
and time averages in k(t) just discussed. A recent Let- 
ter || studies the spectrum of hydrogen in a magnetic 
field. The distribution p{y) is not mentioned, although 
the scaling suggested in Eq.(|^) is tacitly used in taking 
the time average, [over r ~ 0.6]. The results agree quali- 
tatively with the discussion we have presented. There is 
a calculation JhJ of the statistics of the eigenphases of 
an N x N scattering matrix S, for a particular quantum 
chaotic scattering. The function corresponding to k(t) is 

iV -1 |Tr5 tAr | where t is a discrete time, t = n/N, and 
t = 1 is the Heisenberg time. In this case there is an en- 
semble of such systems which was used to help to reduce 
the noise. 

In the limit tw — > oo the noise disappears, and k T (t) 
in that limit becomes self-averaging. For some purposes 
this is adequate. For example Pandey p!l| |, 'proved' that 
'all' spectral correlation functions of particular members 



of an ensemble of random matrices (and thus presumably 
k(t)), are self-averaging [or 'ergodic']. Of course, C(x) as 
defined in Eq. ([!]) is manifestly not self-averaging, since it 
consists of lots of ^-functions, while the ensemble average 
is smooth. But Pandey argued that "observable" quan- 
tities involve an integral over x. Presumably the limits 
of such an integral are fixed and independent of eo. This 
is energy smearing rather than time smearing. Then in 
the limit eo — > oo followed by W — > oo, Pandey indeed 
finds the energy smeared C (x) —* random matrix predic- 
tion. Now one could, if desired, define k(t) as the Fourier 
transform of this C(x). But since the large W limit is 
taken first, for fixed finite x, it's the same as doing a time 
average with tW/A — > oo. 

We are more ambitious. Based on references Q, [[|, 
||, we predict for GUE billiards, that \(k(t)) - k R (t)\ ~ 
for | £ — 1 1 ~ An ensemble average of k is 

indicated and k R is the random matrix form factor. Here 
units are chosen so that A = 1, mass = 1/2, and h = 1. 

We therefore consider the 'renormalized' function 
h(s) = v / ei[fc T (l + s/ v / e^) — k R (l + s/y/eo)] which should 
approach a theoretically predictable limit, 'the signal', of 
order unity, as eo °o provided we can choose a time 
smearing r which leads to a k T (t) sufficiently close to the 
ensemble average prediction. 

Put w ~ eg, and r ~ Cq a ■ Clearly (3 < 1, in order that 
W << eq for large eo . If the time smearing is not to lose 
the signal, then a > \. The 'noise' in this case scales as 
e ( a -/3)/ 2 ^ Multiplying by y/eo, in order to calculate the 
noise in the function h, gives a noise e^ where 7 > ^. 
Thus, we predict that the noise grows at high energy, 
and that the signal will be swamped. 

Perhaps more sophisticated data-analysis schemes 
could extract the signal from the noise, but all those we 
have tried were unsuccessful. One has to keep in mind 
that this is not true noise: for a given spectrum it is 
reproducible. Because, the noise grows at high energy 
in the renormalized h(s), it hurts rather than helps to 
use very high energy data. In favorable cases, it might 
be possible to use low energy data, but that does not 
seem systematically possible. Thus, we characterize this 
situation as giving unexpected numerical difficulties and 
perhaps it is impossible in principle to check the theoret- 
ical predictions, absent an actual ensemble average. 

Fortunately, there are systems which have much larger 
'quasiclassical' corrections to the random matrix results, 
than do strongly chaotic billiards. Most notable is the 
spectrum of nontrivial zeroes of the Riemann zeta func- 
tion. Although the hypothetical 'Hamiltonian' whose 
spectrum coincides with the zeroes is not known, or 
even proved to exist, the 'quasiclassical' parameters are 
known Q. In the language used above, with eo hav- 
ing the mathematical meaning of height in the com- 
plex plane along the critical line, the mean level spacing 
shrinks according to A = 2ir/ ln(eo/27r), which we denote 
A ~ €q°. The 'orbit' periods in units of the Planck time 
are t p ^ m = (mA/27r)lnp, where p is prime and m repe- 
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titions of the fundamental orbit are taken. The stability 
weights, a p in the Gutzwiller formula, are e _mlnp / 2 . We 
expect, based on the techniques of reference |8|, that cor- 
rections of magnitude Cq over a width \t- 1[~ 6q° will 
exist. There should be enough data available to verify 
this prediction. 

Similarly, we think that GUE billiards with strong, 
nonisolated orbits, will also have effects near t = 1 ob- 
servable in principle, but they are probably too difficult 
in practice. The weight of a nonisolated orbit in the 
Gutzwiller formula goes as e rather than e as 
for unstable orbits. We expect corrections of magnitude 
e over a time window e thus allowing r = e ^ 4 . 

— 1/8 

This leads to noise of magnitude e in the rescaled 
h(s). Thus, with many energy levels in the neighborhood 
of the 10 s level, such an effect might be seen in a direct 
calculation of k(t) from the spectrum. In general, a signal 
of width and strength e^" 17 can in principle be extracted 
from the data on a single system only if 77 < g . 

The much studied GOE class of systems (with time 
reversal symmetry) unfortunately has very small effects 
at the Heisenberg time and should be even harder to 
observe. However, GSE systems, (with sympletic sym- 
metry) might be favorable because the RMT result is 
singular at t = 1 and thus could suffer large corrections 
from the nonuniversal effects. These two cases have not 
been worked out in any detail up to now, however. 

Another question is whether it is possible to do en- 
semble averages of the type mentioned, e.g. for a class of 
billiards. This might work if one has an analog computer 
to do the job, namely an experimental system depend- 
ing on a parameter in a nice way, which can be changed 
sufficiently to give a further average. Note that such a 
system is presented in reference ]l0[ . 

We therefore conclude that by manipulations made on 
the spectrum of a specific hard chaos billiard system, it is 
not known how, and it is perhaps impossible, to confirm 
or discredit the theoretical predictions of k(t) for t ~ 1. 
For other systems, particularly the Riemann zeta func- 
tion, it is likely possible. The systems are distinguished 
by the order of magnitude of the nonuniversal effects and 
their dependence on energy. 

Reference [jj) above made no explicit ensemble or time 
averages, but did not find a wide distribution for k(t), but 
rather something like the ensemble average of reference 
Evidently an approximation based on W x /W << 1, 
was made which discarded the noise, so their results are 
correct if the signal has an energy exponent 77 < i. It's 
less clear what reference j| predicts with respect to noise. 
Their result for k(t) is a sum of (5-functions, 5(t — t p ) 
at the orbit periods, but [we argue], these (^-functions 
should have a width A/W. Except for t very small, the 
spacing between the t p 's will be much smaller than this, 
since the spacing decreases exponentially with increas- 
ing t p . [This is the exponential proliferation of orbits in 
chaotic systems.] It was argued 0] that in this case, the 
Hannay-Ozorio de Almeida sum rule can be invoked with 



a smooth, non-noisy result. If that argument is correct 
and carries over to the extensions of the DA used, then 
these methods ||| do not produce the noise. If that is 
so, the noise would have to come from the neglected off- 
diagonal contributions. This might well be the case, as 
it is known that the off-diagonal contributions are very 
important The 'noise' in the Bogomolny-Keating 

formulation thus presents a deep question, which will 
probably not be answered in the near future. 

Certain often used linear transforms of k(t) are self- 
averaging, in agreement with reference pd| ]. These are 
S 2 (L) and A^(L), the number variance and the Dyson- 
Mchta spectral rigidity. The former is given by £ 2 (£) = 
/ dt k(t)[smTrLt/irt} 2 . For L ~ 1 or smaller, there is an 
effective time average over a scale of order unity, thus 
eliminating the noise for large W. For large L, the con- 
tribution of the integral comes from small t, where the 
nonuniversal signal dominates k{t). 
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